GLBRC collaboration¶
Adina, Ashley, Shane, Nejc¶
- Sequencing and alignment Summaries
- Inputs:
- 136 Phred-33 fastq files from samples obtained from washes of Switchgrass and Miscanthus transfered from JGI.
- Illumina TruSeq version 3.2 sequencing adapters.
- Outputs:
- Inputs:
- Overarching analysis - changes in metagenome content over time
- Inputs:
- Functional Gene Count table, mapped to KEGG/functions Annotation Methodology
- Outputs:
- PCoA and analysis of changes over time – related time/season – use
- Synchrony – does the 16S change in concert with functional?
- PROTEST – Procustes superimposition, Mantel
- Assembly of functional genes - is it stable over time, does it accumulate/convergence faster than 16S -John G infant gut - consecutive v. end point –
- Bray-Curtis distance of 16S t1 v t2, T1 v T3, T1 v TFinal, model of convergence for both
- Hierarchical clustering – clusters can inform this question (e.g. T4, 5,6 grouped with final time
- Core functional genes – Nejc is the point person
- Omit core metabolism from the normalized table to focus only on non-housekeeping functions- -
- Discovery/identify Occupancy 100% at a time point
- All genes that are in all samples – Shane is using the median
- Who are those core genes?
- Levy paper – symbiosis genes? How to search?
- Map pfams from the Levy paper to sequence level homology in KEGG?
- Other genes that are not mentioned in Levy but are leaf-specific (e.g., irradiation/desiccation exposure)
- Levy paper – symbiosis genes? How to search?
- Inputs:
- Analysis of bins (MAGs)
- Annotate them – who has which core functional genes (identified above)
- How well do they overlap with our Core from the 16S paper? %coverage – Jackson will be the point person for that
- Match 16S (but this is problematic because probably won't assemble well)
- Use Udeobacter paper as a guide for MAG recovery analysis as compared to 16S – use the most resolved taxonomy –
- Basic stats on MAGs – how abundant are those taxa in the dataset and how persistent are they over time? Where are the gaps, given the core
- metaT
- Map both to metaG and to MAGs
- Subset the core genes for this analysis?
- Do we see the same patterns in metaG and metaT, over time?
- Is everyone active? What proportion of the community? Only certain members? Can we "get away" with just metaG to understand who is active?
Metadata processing for the metagenomic samples¶
#HIDDEN
%matplotlib inline
from glob import glob # nbi:hide_in
from IPython.display import Image
from ipywidgets import Text, Dropdown, Output, interact
from matplotlib.collections import LineCollection
from matplotlib import pyplot as plt
from pandas import DataFrame, merge, read_csv, Series, to_datetime
from pickle import dump, load
from random import randint
import cufflinks as cf
import numpy as np
import qgrid
cf.go_offline()
# cf.set_config_file(offline=False, world_readable=True)
cf.set_config_file(offline=True, world_readable=True)
# For metaG
metadataG = read_csv("metadata/GLBRC_MetaG_Metadata.tsv",sep='\t')
metadataG.set_index("nucleic_acid_name",inplace=True)
metadataG.drop(["source","sampling ID","sequencing_type","height_mean_cm","air_temp_c","rep","Sampling Time","reads_file_name","treatment"],axis=1,inplace=True)
for id in metadataG.index: metadataG.loc[id,"type"] = metadataG[metadataG.index == id].plot_name[0][0:2]
metadataG['Date'] = to_datetime(metadataG.sampling_date)
metadataG=metadataG.sort_values(by=["type","Date"])
metadataG.drop(["sampling_date"],axis=1,inplace=True)
metaG_Reads = read_csv("mapping/metaG/hostRemovalFlagstats/multiqc_data/multiqc_bowtie2.txt",sep="\t")
for id in metaG_Reads.index: metaG_Reads.at[id,"Sample"] = metaG_Reads.at[id,"Sample"].replace(".stat","")
metaG_Reads.set_index("Sample",inplace=True)
metadataG = merge(metadataG,metaG_Reads,left_index=True,right_index=True)
metadataG.sort_values("Date",inplace=True)
metadataG["aligned"] = (metadataG['paired_total']-metadataG['paired_aligned_mate_none_halved'])
metadataG["percPlantAligned"]= metadataG["overall_alignment_rate"]
metadataG["Plant Reads"] = metadataG["aligned"]
metadataG["Total"] = metadataG["paired_total"]
metadataG.drop(['aligned','total_reads','paired_total','paired_aligned_none','paired_total','paired_aligned_one', 'paired_aligned_multi',
'paired_aligned_discord_one', 'paired_aligned_mate_none', 'paired_aligned_mate_one', 'paired_aligned_mate_multi',
'overall_alignment_rate', 'paired_aligned_mate_multi_halved', 'paired_aligned_mate_none_halved', 'paired_aligned_mate_one_halved'],axis=1,inplace=True)
metaG_Reads = read_csv("mapping/metaG/fungalRemovalFlagstats/multiqc_data/multiqc_bowtie2.txt",sep="\t")
for id in metaG_Reads.index: metaG_Reads.at[id,"Sample"] = metaG_Reads.at[id,"Sample"].replace(".fungal.stat","")
metaG_Reads.set_index("Sample",inplace=True)
metadataG = merge(metadataG,metaG_Reads,left_index=True,right_index=True)
metadataG.sort_values("Date",inplace=True)
metadataG["Fungal Reads"] = (metadataG['paired_total']-metadataG['paired_aligned_mate_none_halved'])
metadataG["percFungalAligned"]= metadataG["overall_alignment_rate"]
metadataG['percFungalAligned'] = Series([float("{0:.2f}%".format(val).replace("%","")) for val in metadataG['percFungalAligned']], index = metadataG.index)
metadataG['percPlantAligned'] = Series([float("{0:.2f}%".format(val).replace("%","")) for val in metadataG['percPlantAligned']], index = metadataG.index)
metadataG['Remaining'] = metadataG['paired_aligned_mate_none_halved']
metadataG['percRemaining'] = (metadataG["Remaining"]/metadataG['Total'])*100
metadataG.drop(['total_reads','paired_total','paired_aligned_none','paired_total','paired_aligned_one', 'paired_aligned_multi',
'paired_aligned_discord_one', 'paired_aligned_mate_none', 'paired_aligned_mate_one', 'paired_aligned_mate_multi',
'overall_alignment_rate', 'paired_aligned_mate_multi_halved', 'paired_aligned_mate_none_halved', 'paired_aligned_mate_one_halved'],axis=1,inplace=True)
G5Data = metadataG[metadataG['type'] == 'G5']
G6Data = metadataG[metadataG['type'] == 'G6']
# qgrid_widgetMG = qgrid.show_grid(metadataG)
# qgrid_widgetMG
Metadata processing for the metatranscriptomic samples¶
#HIDDEN
metadataT = read_csv("metadata/GLBRC_MetaT_Metadata.tsv",sep='\t')
metadataT.set_index("nucleic_acid_name",inplace=True)
metadataT.drop(["name", "Year","plotID","month","sampleID","day", "year", "time", "weather", "air_temp_c", "notes", "rep", "time_zone", "longitude", \
"latitude", "altitude", "country", "soil_name", "number_cores", "SPNL_date", "pH", "lime_index", "P_ppm", "K_ppm", \
"Ca_ppm", "Mg_ppm", "organic_matter", "NO3N_ppm", "NH4_ppm", "soil_moisture_percent", "soil_temp_10cm", "LDMC_mg_per_g", \
"nitrogen_percent", "carbon_percent", "carbon_per_nitrogen", "height_mean_cm", "mass_per_leaf_g", "date_of_extraction",\
"nucleic_acid_type", "replicate_extraction", "source", "source_mass", "extraction_method", "elution_vol_ul", "concentration_ng_per_ul",\
"ratio_260_280", "conc_ng_per_g_source", "extracted_by", "sequence_name", "sequencing_date", "conc_sent_ng_per_ul", \
"sequencing_type", "sequencing_facility", "primers", "submitted_for_sequencing", "sequencing_successful", "duplicate_submitted",\
"dup_sequencing_name", "exclude_from_analysis", "itemID_JGI", "sampleID_JGI", "barcode", "JGI_library", "HPCC_path", "JGI_taxonOID",\
"JGI_rawdataname", "time_numeric", "date", "precipitation", "Air_temp_mean", "air_temp_max", "Air_Temp_Min",\
"Air_Pressure","plant_name", "RH", "AH", "Wind_Speed_Mean", "Wind_Direction_Mean", "Solar_Radiation", "PAR", "soil_temp_5_cm_bare_avg", \
"soil_temp_5_cm_sod_avg"],axis=1,inplace=True)
for id in metadataT.index: metadataT.loc[id,"type"] = metadataT[metadataT.index == id].plot_name[0][0:2]
metadataT['Date'] = to_datetime(metadataT.sampling_date)
metadataT = metadataT.sort_values(by=["type","Date","treatment","plot_name"])
metadataT.head()
#Reads After host removal
metaT_Reads = read_csv("mapping/metaT/hostRemovalFlagstats/multiqc_data/multiqc_bowtie2.txt",sep="\t")
for id in metaT_Reads.index: metaT_Reads.at[id,"Sample"] = metaT_Reads.at[id,"Sample"].replace(".stat","")
metaT_Reads.set_index("Sample",inplace=True)
metaT_Reads["percPlantAligned"] = ((metaT_Reads['total_reads'] - metaT_Reads['paired_aligned_none']) / metaT_Reads['total_reads'])*100.0
metadataT = merge(metadataT,metaT_Reads,left_index=True,right_index=True)
metadataT.sort_values("Date",inplace=True)
# metadataT.head()
# qgrid_widgetMT = qgrid.show_grid(metadataG)
# qgrid_widgetMT
Sequencing and alignment Summaries
Adapter Trimming and QC (Trimmomatic) Report
What is the abundance and quality of the reads in each sample?
Step 1. Split JGI files in PE files
split-paired-reads.py --gzip -1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz \$jgi_combined_fastq
Step 2. Use trimmomatic to QC reads
trimmomatic PE -phred33 \$sample.fastq.pe1.gz \$sample.fastq.pe2.gz \$sample.fastq.pe1.gz \$sample.fastq.se1.gz \$sample.fastq.pe2.gz \$sample.fastq.se2.gz \
ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:TRUE \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:36
ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:8:TRUE \
LEADING:3 \
TRAILING:3 \
SLIDINGWINDOW:4:15 \
MINLEN:36
[Outline](#home)
Skip to:
* Adapter Trimming and QC (Trimmomatic) Report
* [Host Plant Alignment Report](#pAlign)
* [Fungal Alignments](#swAlign)
* [Alignment Conclusions](#aConc)
Link to the interactive metaT trim report.
Link to the interactive metaG trim report.
MetaG Trim Report
MetaT Trim Report
Sequencing and alignment Summaries
Host Plant Alignment Report
When the reads are aligned to plant reference genomes, how many of the reads are plant-aligned reads?
Step 1. Align reads to respective host reference assembly
if [[ \$sample == *"G5"* ]]; then
SAMFILE=\$sample.SWGRASS.sam
BAMFILE=\$sample.SWGRASS.bam
bowtie2 --threads \$THREADS -x \$SWITCHGRASS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.SWGRASS.sam 2>\$sample.stat
else
SAMFILE=\$sample.MISCANTS.sam
BAMFILE=\$sample.MISCANTS.bam
bowtie2 --threads \$THREADS -x \$MISCANTHUS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.MISCANTS.sam 2>\$sample.stat
fi
SAMFILE=\$sample.SWGRASS.sam
BAMFILE=\$sample.SWGRASS.bam
bowtie2 --threads \$THREADS -x \$SWITCHGRASS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.SWGRASS.sam 2>\$sample.stat
else
SAMFILE=\$sample.MISCANTS.sam
BAMFILE=\$sample.MISCANTS.bam
bowtie2 --threads \$THREADS -x \$MISCANTHUS
-1 \$sample.fastq.pe1.gz -2 \$sample.fastq.pe2.gz
-U \$sample.fastq.se12.gz -S \$sample.MISCANTS.sam 2>\$sample.stat
fi
[Outline](#home)
Skip to:
* [Adapter Trimming and QC (Trimmomatic) Report](#seqSums)
* Host Plant Alignment Report
* [Fungal Alignments](#fAlign)
* [Alignment Conclusions](#aConc)
#HIDDEN
# %matplotlib widget
from IPython.core.display import HTML
# disp = {"Percentage":{"Switchgrass & Miscanthus": "html/BothPerc.html", "Switchgrass": "html/SwitchPerc.html", "Miscanthus": "html/MiscanPerc.html"}}
# @interact
# def graph(Source=["MetaG","MetaT"],Host_Plant=["Switchgrass & Miscanthus","Switchgrass","Miscanthus"],By=["Percentage","Count"]): #return HTML(disp[By][Host_Plant])
# if Source == "MetaG": metadata = metadataG
# else: metadata = metaT_metadata
# G5Data = metadata[metadata['type'] == 'G5']
# G6Data = metadata[metadata['type'] == 'G6']
# # htmlName = ""; fig = None
# if By == "Percentage":
# if (Host_Plant == "Switchgrass & Miscanthus"):
# # htmlName = "html/BothPerc.html"
# return metadata.groupby(['Date','plot_name']).sum()['percPlantAligned'].unstack().iplot(asFigure=True,title="% Reads Aligning to Respective Plant Host Assembly G5=Switchgrass G6=Miscanthus");
# elif Host_Plant == "Switchgrass":
# # htmlName = "html/SwitchPerc.html"
# return G5Data.groupby(['Date','plot_name']).sum()['percPlantAligned'].unstack().iplot(asFigure=True,title="% Reads Aligning to Switchgrass Assembly",fontsize=12);
# else: # "Miscanthus":
# # htmlName = "html/MiscanPerc.html"
# return G6Data.groupby(['Date','plot_name']).sum()['percPlantAligned'].unstack().iplot(asFigure=True,title="% Reads Aligning to Micanthus Assembly",fontsize=12);
# # ylim = 100; ylab = "% Reads"
# else:
# colors = ["#74C476", "#f26d07","#794955"]
# if (Host_Plant == "Switchgrass & Miscanthus"):
# # htmlName = "html/BothCount.html"
# return metadata.loc[:,['Plant Reads','Fungal Reads', 'Remaining']].iplot(asFigure=True,kind='bar', barmode = 'stack', color=colors, title="Reads Counts By Source")
# elif Host_Plant == "Switchgrass":
# # htmlName = "html/SwitchCount.html"
# return G5Data.loc[:,['Plant Reads','Fungal Reads', 'Remaining']].iplot(asFigure=True,kind='bar', barmode = 'stack', color=colors, title="Switchgrass Reads Counts")
# else:
# # htmlName = "html/MiscanCount.html"
# return G5Data.loc[:,['Plant Reads','Fungal Reads', 'Remaining']].iplot(asFigure=True,kind='bar', barmode = 'stack', color=colors, title="Miscanthus Reads Counts By Source")
# # py.offline.plot(fig,filename=htmlName)
# # return ax #HTML(htmlName)
HTML(filename="html/BothPerc.html")
Sequencing and alignment Summaries
Fungal Alignment Report
What reads are fungal reads?
Step 1. Extract reads that don't align the plant assembly
# A. Convert sam > bam
samtools view -bS \$SAMFILE > \$BAMFILE
# B. Extract all unmapped reads (-f 4) that don't have a mate mapped (-F 256 i.e. both unmapped)
samtools view -b -f 4 -F 256 \$BAMFILE > \$sample.unmapped.bam
# C. Sort the reads
samtools sort -n \$sample.unmapped.bam -o \$sample.unmapped_sorted.bam
# D. Split the unmapped reads in R1 & R2 (back into paired end)
bedtools bamtofastq -i \$sample.unmapped_sorted.bam -fq \$sample.R1.fastq -fq2 \$sample.R2.fastq
samtools view -bS \$SAMFILE > \$BAMFILE
# B. Extract all unmapped reads (-f 4) that don't have a mate mapped (-F 256 i.e. both unmapped)
samtools view -b -f 4 -F 256 \$BAMFILE > \$sample.unmapped.bam
# C. Sort the reads
samtools sort -n \$sample.unmapped.bam -o \$sample.unmapped_sorted.bam
# D. Split the unmapped reads in R1 & R2 (back into paired end)
bedtools bamtofastq -i \$sample.unmapped_sorted.bam -fq \$sample.R1.fastq -fq2 \$sample.R2.fastq
Step 2. Align the reads to the combined fungal assemblies
bowtie2 -x \$FUNGAL -1 \$sample.R1.fastq -2 \$sample.R2.fastq -S \$sample.fungal.sam 2>\$sample.fungal.stat
[Outline](#home)
Skip to:
* [Adapter Trimming and QC (Trimmomatic) Report](#seqSums)
* [Host Plant Alignment Report](#pAlign)
* Fungal Alignments
* [Alignment Conclusions](#aConc)
#HIDDEN
HTML(filename='html/FungalAlign.html')
#HIDDEN
# #HIDDEN
# import plotly.plotly as py
# fig = metadataG.groupby(['Date','plot_name']).sum()['percFungalAligned'].unstack().iplot(yTitle="Percentage of Reads",title="All Samples Perctages of Reads Aligning to Combined Fungal Assemblies",asFigure=True);
# py.offline.plot(fig,filename="html/FungalAlign.html")
# @interact
# def fungalAlign(Samples=["Switchgrass", "Miscanthus", "Switchgrass & Miscanthus"]):
# metadataG.groupby(['Date','plot_name']).sum()['percFungalAligned'].unstack().iplot(yTitle="Percentage of Reads",title="All Samples Perctages of Reads Aligning to Combined Fungal Assemblies")
# if Samples == "Switchgrass": return G5Data.groupby(['Date','plot_name']).sum()['percFungalAligned'].unstack().iplot(yTitle="Percentage of Reads",title="Switchgrass Samples Perctages of Reads Aligning to Combined Fungal Assemblies");
# elif Samples == "Miscanthus": return G6Data.groupby(['Date','plot_name']).sum()['percFungalAligned'].unstack().iplot(yTitle="Percentage of Reads",title="Miscanthus Samples Perctages of Reads Aligning to Combined Fungal Assemblies");
# else: return metadataG.groupby(['Date','plot_name']).sum()['percFungalAligned'].unstack().iplot(yTitle="Percentage of Reads",title="All Samples Perctages of Reads Aligning to Combined Fungal Assemblies");
Sequencing and alignment Summaries
Sequencing Observations
What patterns are there over the season?
- The percentage of plant reads sequenced is high thoughout the season and tapers off towards the end of the season. This may be do to senescence of plant cells.
- Alignment to the 6 selected fungal assemblies is low, but picks up during the warmer months. Overall, reads aligning to the combined fungal assembly were <10% of the non-plant reads.
Annotation Summary¶
The annotations were performed using KEGG's prokaryotic peptides and eukaryotic peptides
Based on the high-overlap, we decided to use the mags instead of the assembly for alignment
MAG Stats¶
checkm lineage_wf -t 20 -x fa Final.contigs.fa.metabat-bins40 metaBinsStats
#HIDDEN
df = read_csv("mags/MAG_Stats.tsv",sep="\t")
df = df.head(100)
df.to_html("html/MAG_Stats.html")
# qgrid_widgetMAGs = qgrid.show_grid(df)
# qgrid_widgetMAGs
HTML(filename="html/MAG_Stats.html")
Overarching analysis - changes in metagenome content over time
All bins with at least 50% completeness were used in a combined assembly. The analysis below comes from the counts for each sample. These counts were normalized by the number of reads that did not align to the fungal or host assemblies.
[Outline](#home)
Skip to:
* Functional annotation accumulation over time
* [Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?](#alphaD)
* [Do we see changes over time in functional richness?](#funcRich)
* [Is there a difference between treatment (Fertilized vs. Unfertilized)?](#treatComp)
* [Is there a difference between host plant?](#plantDif)
* [PCoA of functional counts](#pcoa)
* [Is there similarity between the WGS communities and the 16S data?](#comp16S)
* [How variable are metagenomes across replicate time points?](#permdisp)
COG accumulation over time
Is time a factor that affects the functional richness of phylosphere communities (Alpha Diversity)?
|
|
|
||||||||||||||||||||||||||||||||||||
|
|
|
||||||||||||||||||||||||||||||||||||
Call: adonis(formula = dist.otu ~ map_16S\$time_numeric)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$time_numeric | 1 | 1.9356 | 1.93562 | 47.814 | 0.26298 | 0.001 *** |
| Residuals | 134 | 5.4246 | 0.04048 | 0.73702 | ||
| Total | 135 | 7.3603 | 1.000000 |
Call: adonis(formula = dist.otu ~ map_16S\$treatment)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$treatment | 1 | 0.0289 | 0.028930 | 0.52878 | 0.00393 | 0.755 |
| Residuals | 134 | 7.3313 | 0.054711 | 0.99607 | ||
| Total | 135 | 7.3603 | 1.000000 |
Call: adonis(formula = dist.otu ~ map_16S\$plant)
Permutation: free Number of permutations: 999
Terms added sequentially (first to last)
| Df | SumsOfSqs | MeanSqs | F.Model | R2 | Pr(>F) | |
|---|---|---|---|---|---|---|
| map_16S$plant | 1 | 0.1617 | 0.161704 | 3.0101 | 0.02197 | 0.027 * |
| Residuals | 134 | 7.1986 | 0.053721 | 0.97803 | ||
| Total | 135 | 7.3603 | 1.000000 |
MetaT alignment by count
MetaT alignment by percentage
Next Steps¶
- Fix the metadata to ensure that each sample is accounted for and we have the seq files (Nejc) ✓ Completed Wednesday, Feb 19,2020. Delivered by Kiera
- Modify the R analysis script to read results from the MetaT data (Nejc & Shane)
In order for this to progress, 1 had to be complete. The next step is figure out why the metadata for the MetaT is not meshing with the R script.
- Explore the Kallisto output to determine how to extract annotations (Shane)
Kallisto has an arguement to export a bam file when it runs, but I do not know how these bam files work or if they are the same thing as a bam from bowtie2. I am working with kallisto to get a bam file now and then I will see if that file works to extract annotations from the contigs